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Abstract 

The Residual Smooting Scheme (RSS) have been introduced in [T] as a backward Euler’s 
method with a simplified implicit part for the solution of parabolic problems. RSS have 
stability properties comparable to those of semi-implicit schemes while giving possibilities 
for reducing the computational cost. A similar approach was introduced independently in 
[101 111] but from the Fourier point of view. We present here a unified framework for these 
schemes and propose practical implementations and extensions of the RSS schemes for the 
long time simulation of nonlinear parabolic problems when discretized by using high order 
finite differences compact schemes. Stability results are presented in the linear and the 
nonlinear case. Numerical simulations of 2D incompressible Navier-Stokes equations are 
given for illustrating the robustness of the method. 
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1 Introduction 

It is a common fact in numerical analysis that the choice of a time marching scheme must 
balance stability, accuracy and reasonable computational cost. Typically, when considering e.g. 
the numerical solution of a space-discretized parabolic problems, such as 

^+Au = f, 
u{0) = uo, 

where A denotes the stiffness matrix, it is well known that the implicit schemes are stable but 
need an additional problem to be solved at each step while the explicit schemes are very cheap 


1 


but suffer of a hard time step limitation making them bad suited for capturing the long time 
behavior of the solutions. An ideal scheme should combine stability and low computational cost 
(explicity) for a comparable accuracy. 

Two independent attempts have been made successfully in that direction : 

First, A. Cohen and al proposed in [1] the following stabilization to the forward Euler scheme 
(Residual Smoothing Scheme or RSS) for the discretized parabolic equations associated to ho¬ 
mogeneous Dirichlet (or Neumann) Boundary conditions: 

'- V -^ (^j 

Stabilization term 

where r is a positive real number to be chosen and B a preconditioner of the stiffness matrix 
A. Originally introduced in the context of wavelet discretizations, the matrix B can be taken as 
the diagonal part of A and is then an inconditional preconditioner: the new scheme is no more 
expensive than the classical forward Euler’s while the stability is increased. However, RSS is 
only first order accurate in time and, in order to increase the accuracy, it was proposed in [T] to 
apply a Richardson extrapolation; typically a second order of accuracy was obtained as shown 
by numerical evidences. A rough analysis of the stabilized and extrapolated scheme was made 
by Ribot and Schatzman [iiiin!, they derived stability and error estimates in energy norms . 

Independently, Costa, Dettori, Gottlieb and Temam have introduced in [101 m a similar 
approach but starting from a Fourier-analysis point of view, in the context of multiresolution 
methods of nonlinear Galerkin type for spectral discretizations (Fourier, Chebyshev), see |19j . 
They proposed to stabilize the forward Euler scheme for the heat equation 


Uxx — f : 

by adding a stabilization term of the form The new scheme writes in Fourier basis as 


u 




= + /m, rn = Nj_i, • • • , Nj, 


Stabilization term 


(3) 


where we have decomposed the frequency range [1,N] into A(j]. In other words 


^{k+l) _ ^{k) 


At 


- R(#+^) - + /, 


with D = diag{l,4:, 9, • • • and 


B = 


( Bi 0 
0 

Vo 


0 \ 


Bd J 
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with Bj = I3jldj, Idj being the identity matrix of size {Nj — Nj-i) x (Nj — Nj-i) (we have set 
Nq = 0 and = N)-, B is then a preconditioner of D in the Fourier space. In the linear case, 
these two approaches coincide. Of course, same framework can be derived when considering 
orthogonal polynomials. Costa and Chehab mi have extended this scheme to hierarchical 
discretizations in finite differences. 

The main advantage of the RSS approach is that a simplified (yet costless) solver is used for 
the implicit part of the time marching scheme while displaying comparable stability properties 
to Backward’s Euler Scheme. One situation of particular interest, on which we focus in the 
present work, occurs when handling high order discretizations of the stiffness matrix A, e.g. 
with finite differences compact schemes. In that case, A is full, this is due to the implicit part of 
the scheme. Hence, matrix-vector product are costlty and must be reduced as far as possible. A 
lower level of space dicretization (say second order) generates a sparse stiffness matrix B which 
is a natural efficient preconditioner of A. Then, the RSS scheme can be implemented efficiently 
taking advantage of the existing (sometimes fast) solvers of the system of the form 

{Id + AtB)u = /, 


such as sparse factorizations and FFT. 

The RSS approach can be proposed to solve fully discretized time dependent PDEs with high 
accurate spatial discretization, with compact scheme, while using the computational facilities of 
the sparse numerical solvers (Fast solvers, limited memory). 

In this article, we propose a unified approach to RRS-like schemes that rely [I] and m- We 
derive stability results in the linear and the nonlinear case; we also present practical and efficient 
adaptations for the high accurate finite differences solutions of nonlinear parabolic equations. 

The paper is organized as follows: in Section 2, we derive a general approach to RSS schemes 
and we give stability results in the linear and the nonlinear case, the accuracy in time of the new 
scheme is also discuted. After that, in section 3, we describe the compact scheme discretization 
and the preconditioning that will be used. Section 4 is devoted to numerical illustrations, 
we compare RSS approach to the classical one with emphasis on the stability, the accuracy 
(particularly the dynamics to steady state) and for that purspose, we solve high accurate finite 
difference steady state of 2D incompressible Navier-Stokes equations (Lid driven cavity) and 
recover the results of the litterature. 

2 Derivation of the stabilized schemes - properties 

2.1 Formal Derivation of the stabilized schemes 

Let us consider the finite dimensional differential system 

du ^, 

— +F{u) = 0,t > 0, 
u{0) = uo, 
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(4) 

(5) 


here F : IR^ —)• ]R^ is a regular map. The backward Euler scheme applied to the above system 
generates the iterations 

y{k+i) _ ^(k) ^ = 0, 

and a (possibly) nonlinear problem must be solved at each step. Making the approximation 

where F'{u^^^) denotes the differential of F at we obtain the scheme 

^ - uW) + = 0, 

so 

^(k+1) ^ ^(k) _ At(/d + 

Setting <l>(u) = v — + AtF(v), we see that is nothing else but the first iteration of the 

Newton-Raphson scheme applied to $(u) when starting from the initial guess 
Now, if we replace by a preconditioner tB^, we find 

^(fc+l^_^(fc) + +F(uW)=0, 

Global stabilization 

and is thus the first iteration of a quasi Newton Method applied to d>(u) when starting 

from the initial guess 


The efficiency of this stabilized scheme is closely related to the cost of the computation of 
the preconditioner of the Jacobian matrix which changes at each iteration: technique of existing 
updating factorizations as those presented in [6] and [2] could be adapted. 

In the present work, we will not discuss on the analysis of the nonlinear version of the scheme, 
say but we will present on Section 4 numerical results obtained with this scheme. We will 
rather consider the semi linear approach: if F{u) can be expressed as F(u) = Au + f{u), we 
define teh scheme 

^(fc+i) u(k) ^ _y(fc)) 

'-V-" (7) 

Stabilization of the linear part 


where R is a preconditioner of A 

It is important to point out that 0 is consistant with the computation of steady states and 
can then be applied as a pseudo-time numerical solver, as illustrated in Section 4 with impress¬ 
ible NSE. 
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Of course this stabilization approach applies to the linear case (/(«) = 0). Particularly, RSS 
is a simplified 0-scheme in which the matrix A is replaced by a preconditioner. Indeed, the 
0-scheme write, after simplifications as 

Jk+i) ^ ^(k+i) _ ^ eAtA)-\Au^’^'> - f) 

and, substituting A by R in the implicit part, we recover the RSS 

^(fc+i) ^ ^(fc+i) _ + eAtB)-^{Au^'^^ - /) 


with 9 = T. 


In practice, the use of a preconditioner B oi A leads to propose K = Id + rAtB as pre¬ 
conditioner of M = Id + AtA, where Id is the identity matrix. This can be realized in many 
ways, e.g., by computing K as an incomplete factorization of M; in some cases it can be done 
by solving the linear systems involving K with fast solvers (FFT or so), see section 3. The RSS 
approach applies also to linear problems with a matrix A{t) which depends on time t: 


du 

dt 


+ A{t)u = /, 


that we discretize as 


(Id + rAiBk) = F — A{kAt)u^^'>' 


K. 


( 8 ) 


(9) 


The matrix can be computed as an incomplete LU factorization of M = Id + AtA{kAt) 
and, if A{t) does vary slightly with t, incremental factorization updates from K^-i can be done 
following the techniques proposed in [6]. Notice also that scheme (10) can be obtained by 
applying RSS to linearized equation, as 


{Id + rAtBk) 

'-V-" 

Kk 


( 10 ) 


where Bk is here such that F{u^) = BkU^ 


2.2 Properties of the schemes 
2.2.1 The linear case 

Let A and B be both N x N real symmetric positive definite matrices; the symmetry of A is 
considered for the sake of simplicity however the following approach remains valuable in the 
nonsymmetric case, see section 3 and Theorem 3.2, . We assume that there exist two strictly 
positive real numbers a and j3 such that 


n 


a < Bu,u ><< Au,u >< j3 < Bu,u >, Vtt G 
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It is important to note that a and [3 can depend on the dimension N, if not the matrix B is said 
to be an inconditional preconditioner of A. We will use the following notations: < . > is the 

euclidian scalar product in and || . ||, the associated norm. We will note Xmin (resp. Xmax) 
the lowest (resp. the largest) eigenvalue of A. 


We now consider the RSS scheme applied to the discretized heat equation 

„,(A:+ 1 ) (fc) 

^ - uW) = 

At 

We first prove a simple stability result: 

Proposition 2.1 Under hypothesis T-L, we have the following stability eonditions: 

• IfT>^, the scheme is unconditionally stable (i.e. stable V At > Oj 

0 2 

• If T < ^, then the scheme is stable for 0 < At < -? --. 


Proof. Taking the usual scalar product of each terms with we find 


1 

At 


fk+i) _ ^(k) ||2 ^ >= - < > . 


Using the parallelogram identity, 

< >= ^ (< >-< A(u(^\u(^^+ < A(u^^+^^ > 

we infer 

^ II f +r < > 

-2 (< >) = 0 . 

Hence the stability condition < > holds when 

^ II f +T < » 0. 

We have, using % 

T < > “2 < > 

> 

— 2 ^ < > (> 0 ). 

A sufficient stability condition is then 

^ II f + (^ " 0 < - nW » 0. 
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This is satisfied once + Xmin{^ “ 2^ ^ Therefore, if ^ j ^ 0 the previous inequaiity 
hoids for aii At > 0, this means the stabiiity VAf > 0. 

Now if r < ^, then, since < >< p{A) || |p, a 

sufficient condition of stabiiity is 

+ (^ - 1 ) p(^) > 0. 

from which we deduce 

2 

0 < At < 7-r-, 

(l-2j)p(A) 

as a sufficient stabiiity condition. ■ 

We point out that if i3 = ^ (then a = /3 = 1) and t = 9 € [0,1], the stabiiity condition coincide 
with the one of the 0-scheme. 

The stabiiity is easiiy obtained when taking r iarge enough. However, a too iarge vaiue of r 
deteriorates the consistency of the scheme and, as a particuiar effect, the convergence to the 
steady state is ionger in time. In fact both the vaiue of r and the preconditioning quaiity of 
B act on the accuracy of the RRS scheme which remains aiways first order accurate in time as 
iiiustrated in Figure ([^. 




Figure 1: (ieft) Influence of the parameter r on the accuracy of RSS. Error vs time for different 
values of r - iV = 127, At = 0.004, (right) Error vs time with ILU(2) and SSOR preconditioners, 
N = 63,rl. At = 0.004 


A first natural question is the choice of the best value Topt of r, for a fixed preconditioner B] 
Topt can be simply computed such as minimizing || tB — A Hi?. We easily show that 


'^opt — 


trace{BA' A) 


B 


|2 

If 
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We remark that when B = Id, then Topt = the mean value of the eigenvalues of A. 

A second natural question deals with the gain of stability brought by the RSS scheme as 
respect to an explicit method, namely forward Euler’s for which the time step must be taken 


strictly lower than Ate = 


P(A)- 


In other words, for a given preconditioner matrix B and for a 


given number r, we look to k > 1 such that RSS is stable with time step At = HiAtc- We deduce 
directly from the previous computations. 


Proposition 2.2 


R 

If T > RSS is infinitely more stable than backward Euler’s. 




then RSS is at least k = 


1 


Proof. We deduce from the proposition 


l-W 

2.1 that 


times more stable than Euler’s. 


1 - ^ ) Pi^) 


= K 


P{A) 


,then K = 


1-^ 




We now propose to quantify the consistency error with r by comparison with the Backward 
Euler Scheme (which is unconditionally stable). Particularly we analyse the behavior of the 
difference of the sequences generated by the two schemes: the stabilization has as an effect to 
slow down the convergence in time to the steady state. 

Proposition 2.3 We consider the two sequences 

r. , (/C + 1) , (/C) 

^ ~ ^ - nW) = / - 


and 


with We let M = 

there exists 7 G [0,1[ such that 




At 




Id — At{Id + rAtB) ^A 


= /, 


and we assume that 


M II< 1, then, 


u(k) _y(k) II< ^^2 \\yB-A 


1 


1-7 


/ — II, VA: > 0. 


Proof. We first remark that || M ||< 1 implies the stability of the RSS scheme, since M is the 
iteration matrix and p{M) <|| M ||. 

We take the difference and we let . We have 

^ + {tB - A){v^’^+^^ - = -Aw^^^ 

Hence, after the usual simplifications, we can write 

yjik+i) ^ ^ TAtB)-^A)w^'^^ - At{Id + TAtB)-\TB - A){v^’^+^'> - 














Using the definition of M, we obtain directly the estimate 

II ||<|| M III! || +At || (Id + rAtB)-^ |||| tB-A\\\\ t;U+i) _ 

We first remark that we have the relations 

y(k+i) _ yik) ^ ^ AtA)-^{f - Av^’^'>), 

and 

^(k+i) ^ _ ^tA){Id + AtA)-^A’^\ 

where we have set = f — Av^^\ It follows that 

^(fc+i) _ y{k) ^ ^ AtA)-^ {{Id - AtA){Id + AtA)-^ fA°\ 

Then 


II yik+i) _ yik) ||< II ^ AtA)-^ IIII {Id - AtA){Id + AtA)-^ f II || . 

We set 7 =|| {Id — AtA){Id + AtA)~^ ||, we have of course 7 < 1, VAt > 0. A simple induction 
gives 


k 

II ||<|| M f II || +AA || {Id+rAtB)-^ |||| tB-A |||| {Id+AtA)-^ || E 7(-?') II M \\^~^ II r® 

i=o 

Using II A/ ||< 1 and = 0 we find 


k 

II yjik+i) II < II {Id + rAtB)-^ III! tB-A\\ At || {Id + AtA)-^ || ^7^^'^ || || . 

j=0 

Hence 

II ||< AA II tH - a II II II, 

1-7 

which shows the dependence on || tB — A ||. Of course if r = 1 and B = A we have = 0,Vfe. 
Hence the result. ■ 


We deduce immediately the 

Corollary 2.4 The RSS method is first order accurate. 

Proof. It suffices to own that the Backward Euler method is first order accurate. ■ 


We now will consider the particular the case in which H is a diagonal matrix. This choice 
allows a fast solution of the implicit part of RSS, the matrix Id + rAtB being also diagonal. 
In general situations, a diagonal preconditioner is not the most efficient, but in some cases, e.g. 
when the discrete problem in written in hierarchical-like bases, it is particularly interesting to 
consider B as a diagonal matrix: A. Cohen et al [T] have introduced the RSS scheme for a 
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problem discretized in a wavelet basis, in which the diagonal part of the stiffness matrix is an 
inconditional preconditioner; this was independently applied by Costa et al san] using Fourier 
and Chebyshev expansions and then in finite differences with incrementals unkowns by Chehab 
and Costa laisii- The underlying idea is to decompose the unknowns U of the nodal basis 
into a hierarchy of arrays of details at different levels {Uo, C/i, • • • , here Uq is associated to 

a coarse discretization and then captures only low frequencies while the other set of components 
Uj, j = 1, • • • ,d are details associated to refinment of the coarse approximation space and cap¬ 
ture high frequencies. 

The time limitation of the explicit schemes for parabolic problems depends on its capability 
to contain high frequencies expansions. In Fourier-like basis, the details attached to high fre¬ 
quencies are small quantities regardless to the details attached to low frequencies since they 
contribute residually to the energy norm of the signal. As proposed in [miaiE] this situation 
allows to damp differently the high and the low frequencies components without deterioring the 
consistency of the scheme. 


Consider for instance the heat equation 

Ut Uxx — f 1 

that we discretize in time with the forward Euler scheme as 

(A:-1-1) _ (fc) 

_ _ _ _ ,,(fc) _ 


At 


= /■ 


- = /• 


In [TDl m] , this scheme was proposed to be stabilized as 

-+ ^- 

Considering a Fourier discretization, we find, after simplifications 

^(k+i) 4 k) o 

+ fm, m = 1, • • • , A. 

If we decompose the frequency range [1, N] into Nj], the above scheme can be applied 

to each range with a different stabilizing parameter /3j, so we obtain, for j = 1, ■ ■ ■ ,d 

4k+l) _ 4k) n 

+fm, m = Nj.i, ■ ■ ■ ,Nj. 

Letting Pi = the last scheme is rewritten as 


u, 


(fc-i-i) — / 1 _ 


m?At 


-L 


At 


1 -|- AtPi J 1 -|- Atp,, 

The stability condition is then 

{m? — 2pj)At < 2,m = Aj-i, • • • , Nj 


fm, rn = ,Nj. 
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Ni 


The scheme is unconditionally stable at level j if and stable under condition 0 < 

2 

At < — 5 ^ otherwise. This condition is of course similar to that found in Proposition 


EH 


iV/+i - 2/3, 


Now, considering all the components, we write 

- + f, 


uik+i) _ 


At 

with D = diag{1^4:, 9, • • • N"^) and 


B = 


( Bi 
0 

Vo 


0 \ 


Bd J 


where Bj = I3jldj, Idj being the identity matrix of size {Nj — iV,_i) x (Nj — lV,-i) (we have set 
A^o = 0 and Nd = N). 

B is then a preconditioner of D in the Fourier space; in the linear case, these two approaches 
coincide. Costa and Chehab [Hi have extended this scheme to hierarchical discretizations in 
finite differences. 

Note that this framework allows to damp high frequencies and to leave unchanged the low ones 
changing only slightly the consistency of the scheme while increasing its stability. This can be 
done typically by taking jdj = 0 on the low freqencies components and /?,• >> 1 on the high ones. 

As stated in the introduction, we concentrate on problems discretized in the nodal basis, 
however, it is useful to make the link with the hierarchical approach. We give a block version 


of Proposition 2.1 


We assume that the stiffness matrix A written in a detail basis posseses the following block 
decomposition 

/ ^ 1,1 Ai^ 2 ■ ■ ■ Ad,d \ 

^ 2,1 

A = 

\ • • • Ad^d j 

We note the corresponding block decomposition of a vector U as U = (ui, • • • ,Ud)'^- We have 
the 

Theorem 2.5 A sufficient stability condition is 


^ + T/3i - - II Aij II + II A,- i II j 


> 0, i = 1, • • • ,d. 


Proof. Taking the scalar product of the equation with we hnd 

^ II - C/W f +r < > 


-^ (< AU^^\ >) = 0 , 
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and using the block decomposition 


d 

E l II (fc+l) (k) ||2 I n II (fc+1) (k) ||2 

— II ’ -u\’ II +r/3i II u\ ’ -u\> II 


i=l 


d d 


(< u^k+i) 


i=l j=l 


A sufficient condition for the energy stability is then 


E 


^ II Il2 


At 


yi^) iii! I a II (*^+l) {^) ii2 ^ 

- u\ Ip +r/?i II u\ -u\ ’ Ip-- 


d d 




2=1 
Since, 

1 


i=l j=l 


d d 


d d 


2^'^^ AiK' 

i=l j=l 


(fe+1) (fc), (fc+1) (k) 

— Uj ),ul — u 


f > £- 1 EEM‘. 


J III 


(fc+1) (fc) II (fc+1) (fc) 
- Uj- 11 u • -u\ ' 


i=l 7=1 

d d 


s-lEEM.. 

2=1 ?' = 1 
d d 

II 

7=1 j = l 


_(fc+l) ^ (fc) |2 


'-Uj 


u. 


(fc+l) ^(k) |2 


— U} 


We find as sufficient stability condition 

d 


s(i 


+ J I X/ II II II "^*’1 


'12^ Y1 II ^*-1 II + II 

7=1 \i = l 


(fc+l) (fc) ||2 ^ ^ 
111 — u) > 0 . 


vi=i 


In particular, if Tj3i > ^ || Aij || + || Aij ||^ >0, the stability is unconditional. ■ 

It has to be noted that when d = 1 and B = Id, we recover the result given by Proposition |2.1[ 

We infer also that, the extradiagonal coefficients of stiffness matrices in hierarchical-like 
basis enjoy of a descreasing magnitude properties far from the diagonal making successful the 
approach. 

If we consider Fourier basis, the stifness matrix A is diagonal and the stability condition at 
range j is 

(m^ — 2Pj)At <2,m = Nj^i, ■ ■ ■ ,Nj 

iV 2 1 

so, taking 13j = —we have an incondiditionally stable RSS scheme. 
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2.2.2 The nonlinear case 


We aim at applying the RSS Scheme to reaction-diffusion equation (such as Allen-Cahn’s, see 
section 4) say 


Art -|—2 ~ X G n, t > 0, 

(11) 

1^=0 9!2,t>, 

on 

(12) 

u{x, 0) = xxo(x) X G II. 

(13) 


The RSS scheme applied to the discretized scheme reads as 
„,(fc+l) _ „,{k) 1 

- -= -AnW (14) 

We set E{u) = ^ < Art, u > < F{u), 1 >, where F is a primitive of / that we choose such 

that F(0) = 0. We say that the scheme is energy descreasing if 


Particularly, if F as nonnegative values, the scheme will be stable for the norm | u \a= 
^/< Au, u >. 

Theorem 2.6 Assume that f is and \ f |oo< L- We have the following stability conditions 
• If T > ^ then 


— “ 2 j — 2 ^ — 0 then the scheme is unconditionally stable 

— “ 2 ) “ 2 ^ ^ then the scheme is table for 

1 


L 


0 < At < 

If T < ^ then the scheme is table for 

0 < At < 






^ - (J - 2) 


Proof. We have, taking the scalar product in IR^ of each element of the with 

^ II f +T < R(u(^+^) - >= -^ < > 

— < — xx^^) > . 
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Let us now consider the nonlinear term. We have using Taylor-Lagrange’s expansion 




Taking the sum of these term for i = 1, • • • ,N, we obtain 

so that 


N 


i=l 


>< -^(< >) + ^ II f . 

The other terms are treated exactly as in the linear case, and we use the identity 

< >= ^ (< > 
Making use of these results, we find, after the usual simplifications 

(it " f ^ > +E{u^’^+^^)-E{u^^^) < 0 

Hence the stability is obtained when 

(it " ^ (i " 0 ^ » 0 . 

Hence the result. ■ 

We point ou that in the linear case (/ = 0 and L = 0), we recover the stability conditions 
given by poposition |2.1[ 

As an application, we can consider the simualtion of the Allen-Calm equation 

du 

m 


An - 1 - \f{u) = 0 

X £ Q,t > 0, 

(15) 

du 

dn 

X £ dQ,'it > 0, 

(16) 

u{x, 0) = uo{x) 

X £ Q. 

(17) 


This reaction-diffusion equation describes the process of phase separation in many situations, 
[?]. In practice, we will chose f{u) = uiv? — 1). 

Notice that Shen et al [l8] have proposed the scheme 


„,(fc+l) _ „,(fc) C* 1 

- = 0 . 


(18) 


With Theorem 2.6, we recover the stability conditions proposed by J. Shen when A = B and 
r = 1 and 5 = 0; The term plays the role of the stabilizator, following the 
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same principle as the schemes introduced in [lOL ITT] . The time restriction become harder when 
e takes small values. This situation motivates the use of stabilized schemes. Now, if S' > the 


scheme (18) is unconditionally stable. This is to be compared with the RSS scheme (14). We 


find as inconditional stability condition 


r > 


^ f3e 


which is a comparable condition for e small enough and /3 bounded, since in pratice Xmin is a 
positive constant which not depend on the dimension of the problem, the first eigenvalue of the 
stifness matrix is indeed nicely captured by the discretization schemes. However the additional 
stabilizing terms can deteriorate the consistency. 


2.3 Richardson extrapolation 

In [T] the authors have proposed to increase the accuracy of the stabilized scheme by smoothing 
the residual using a Richardson extrapolation process: 

The solution of 


du 

dt 


F{u), 


by the forward Euler scheme dehnes the iterations 


^k+i = AtF{u^) = 


The smoothed sequence is defined by 


VI = GAtiu'^), 

V2,0 = GAt/2{u^), 

V2,l = GAt/2iv2,o), 

y^k+1 _ 2 v2^i — Vi- 

It is second order accurate in time. The accuracy of the stabilized scheme is increased by 
applying the extrapolation to the iteration operator 

GAtAu^) = u’‘ + At{Id + TAtB)-^F{A). 

The improvement is studied analytically in IIZI, but numerical evidences point out the efficiency 
of the approach, see also the numerical results presented in Section 4. 

Below the Extrapoled RSS scheme 


15 





Algorithm 1 : Extrapoled RSS Scheme 

1 : given 

2 : 

3: for k = 0,1, ■ ■ ■ until convergence do 
4 : Solve {Id + TfB)vi = -fF{u^^) 

5: Set ui = + vi 

6 : Solve {Id + t^B)v 2 = —^F{ui) 

7: Set U 2 = Ui + V 2 

8 : Solve {Id + TAtB)v 3 = —AtF{u^^^) 

9: Set U 3 = + V 3 

10 : Set = 2 u 2 — U 3 

11 : end for 


3 Discretisation in space and preconditioning 

3.1 Compact FD Scheme Discretization 

A way to obtain a high level of accuracy with a finite difference scheme, that can be compared 
with the spectral one, is to implement finite difference compact schemes M- These schemes 
consist in approaching a linear operator (differentiation, interpolation) by a rational (instead of 
polynomial-like) finite differences scheme. We describre briefly here only their construction for 
the approximation of the first and the second derivative and we refer to m for more details. 
We first consider the schemes in space dimension one. 

Let U = {Ui, - ■ ■ ,Un)^ denotes a vector whose the components are the approximations of 
a regular function u at (regularly spaced) grid points Xi = ih, i = I,-- - ,n. We compute 
approximations of 14 = C{u){xi) as solution of a system 

P.V = QU, 

so the approximation matrix is formally B = P~^Q. When P = Id, the scheme is explicit 
and we recover the framework of classical finite difference schemes; when P ^ Id, the scheme 
is implicit an dit is possible to reach high order accuracy, the implicity allows to mimmic the 
spectral global dependence. In practice, P is a banded sparse matrix easy to invert and very 
well conditioned making the compact scheme approach robust and not costly regardless to the 
precision brought. We here give the matrices P and Q for the fourth order approximation of 
the first and the second derivative, details can be found in |14j . 

Let us begin with the first derivative. We have 

/I i \ 

i 1 

p — 4 

V 4 h 
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/ ai 


2h 


V 


02 

0 


— 04 


03 

3 

2 


_3 

2 

-03 


04 


0 

-02 


-Ol/ 


with Ol = —2, 02 = 3, 03 = — i and 04 = ^. 


In the same way, we can build the fourth order compact schemes for the second order derivative. 
We first consider the compact scheme asociated to the discretization of the second derivative 
with homogeneous Dirichlet boundary conditions. We have 
with 


/I 

10 


P = 


10 

1 


and 


/ Ol 




02 

12 

5 


03 

_6 

12 ^ 


10 


10 


O4 


1 

_!_ 

10 


\ 


10 

1 / 


V 


OAr-4 CLN-S 


O5 


12 

5 

07V-2 


12 

5 


ttN-l aN / 


here the constant oi, 02 , 03 , ... are given by 


Ol 

02 

03 

04 

05 


Oat 

OAT-l 

OAr-2 

OAr-3 

OAr-4 


_67 

(^ 0 ’ 

“T 2 ’ 

13 

10 ’ 

61 

^ 120 ’ 
12 - 


Now applying the same approach, we can consider fourth order compact schemes for the 
second derivative with associated homogeneous Neumann Boundary conditions 
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and 


with 


^ h 2 


ai 

_6 

5 


02 

12 

5 


03 

_6 

ll 

5 


04 

_6 

5 


O 5 


_6 12 _6 

5 5 5 

_6 12 

5 5 

OAr-4 OAr_3 a]Sf-2 O 7 V-I 


_6 

5 

otv/ 


oi = (In 
02 = a-N-i 

< 03 = aAr_2 

04 = aN -3 

, 05 = aAr_4 

To obtain the compact schemes of first and second order derivative in space dimension 2 and 
3, it suffices to use the previous schemes and to expand them tensorally. 

The finite differences discretization matrices of derivative on cartesian domains can be obtained 
by those on the interval using Kronecker products. Indeed, if denotes the discretization 
matrix on [0,1] associated to Dirichlet Boundary conditions, using N internal discretization 
points, then 


2681 

^ 1 ’ 

111 ’ 
40 ’ 
_13 
15’ 
59 

480- 


/dM (8) 


will be the discretization matrix of the same operator but on a grid composed oi N x M points, 
the corresponding laplacian matrix will be A^^®IdN + IdM®A^^. We recall that the Kronecker 
product of a m X n matrix Ahy a, p x q matix B is defined as 


A^B 


( aiiB 
\(lmlB 



3.2 Preconditioning FD Compact schemes using second order FD 

The matrices associated to compact finite difference schemes are full, this is due to the implicit 
nature of the scheme. A natural idea to built a sparse preconditioner is to use the matrix obtained 
by applying a lower accurate finite discretization scheme, particularly a second order one. We 
here describe the approach for one dimensional problem, extensions to higher dimensions are 
obtained using kronecker products, also we restrict to linear problem for simplicity. Let A 2 (resp 
A 4 ) be the second order (resp. the fourth order) discretization matrice of — A on a regular grid 
composed of N internal points. The RSS scheme writes as 

~ + A 4 u(^) = /. (19) 
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The numerical treatment of non homogeneous (possibly time depending) Dirichlet boundary 
conditions can be realized with the RSS approach. Indeed, let us note Am{u,n), m = 2,4, the 
mth order finite difference discretization of —A of u with Dirichlet conditions at time nAt, note 
that this operator is affine. The stabilized scheme writes formally as 

(A:+l) (fc) 

- -+ T{A2{u^'^+^\k + 1) - A2{u^^\k)) + A^{u^^\ k) = /, (20) 

Making the approximation A 2 {u^^^^\ fe + 1) ~ A 2 {u^^~^^\ k), we obtain 

- u(^)) + A^iu^^lk) = /. (21) 

It is to be pointed out that for the solutions of 2D and 3D Poisson problems, the number of 
iterations to the convergence is not dependent on the dimension of the system. Also, in the case 
of Poisson-like problem, we can use FFT to solve the preconditioning system, speeding up the 
resolution of the original linear system. This approach will be followed also for the fast solution 
of the heat equation that arises in parabolic problems. 

Remark 3.1 A analogous approach have been done in the context of hierarchical preconditioners 
in finite differences, where the fourth order discretization matrix of —A was preconditioned by 
the hierarchical transfert matrix attached to the second order accurate discretization of —A, see 

m- 


Below, we report numerical results on the solution of 2D and 3D Poisson problems when dis¬ 
cretized by fourth order compact schemes. The preconditioning matrix is obtained by applying 
corresponding second order finite difference schemes. We took a random r.h.s b=l-2=i'rand(N, 1) 
so that many frequencies including high ones are present. The initial data is u = 0, the tol¬ 
erance parameter 10“^^. The result we report is the maximum number of external iterations 
to convergence, on 5 independent numerical resolutions, the number of discretization point per 
direction n is reported as (n). The stiffnes matrices are then of respective sizes x (2D 
problem) and x (3D problem). 


Problem 

# it- (n) 

if it. (n) ) 

if it. (n) 

if it. (n) 

#it. (n) 

#it. (n) 

Poisson 2D 

12 (n=15) 

11 (n=31) 

10 (n=63) 

10 (n=127) 

9 (n=255) 

8 (n=511) 

Poisson 3D 

12 (n=15) 

11 (n=31) 

11 (n=63) 





Table 1: Solutions of 2D and 3D Poisson problem with GMRES and second order preconditioner 


Here, the fourth order discretization matrix A of —A is nonsymmetric while the precondi¬ 
tioning matrix B is. However, in practice the RSS method is still efficient. This is due to the 
small size of the skewsymmetric part of A. Indeed, denoting by <5 =|| A — ||, we can prove 

the following general stability result. 

Theorem 3.2 Let A G A4,i(M'^). IFe assume that A is positive definite and B a symmetric 
definite positive preconditioning matrix of A satisfy hypothesis %. We set <5 =|| A — A'^ || and 
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= (/3^ — 2ar)f + ^(5^. Assume that §■ --— , ^ > 0. Then the RSS scheme has the 

following stability conditions 
/3^ 

i- if T> T^ + -—then the scheme is unconditionally stable. 
oCxAminK-^ ) 


ii- Ifr < §-- 


mm \ 
6^ 


2 ck 8a\max{BY 


then the scheme is table under condition 


0 < At < 


2a 


•^(XmaxiB)) 


Q'2. r2 

Hi. If Tj- — -—< T < if- + - —r——TT^ then the scheme is table under condition 
^ 8aXmax{B)^ - 8aXmin{BY 


/3" , 52 


0 < At < 


2a 


<I>{Xmin{B)) 


Tf 


<52 


<T < iYr - 


6^ 


2a 8aX^in{BY ^ ^ 2a 8aX^ax{B? 


then the scheme is table under condition 


0 < At < 


2a 


Max($(Amm(-6))) ^(A max (B))) 


Here Xmin{B) (resp. Xmax{B) denotes the lowest (resp. the largest) eigenvalue of B. 

Proof. The RSS scheme reads as 

^(fc+i) ^ ^(k) _ ^ T6tB)-^ 

and is stable under the necessary and sufficient condition p{M) < 1; in the general case the 
eigenvalues of M can be complex. Let v G be an eigenvector of M associated to the 
eigenvalue A = a + ib. We have 

Xv = Mv, so (1 — A) {Id + rStB) v = AtAv. 


We decompose u as u = ui + iv 2 and, separating real and imaginary parts, we obtain 
Avi = (ui + rAtEvi) + {^2 + TAtBv 2 ), 


Avi = 


1 — g 


We have then, on the one hand 


< Avi,Vi > + < Av2,V2 > = 
and, on the other hand. 


1 — g 
At 


¥ 

{v 2 + TAtBv 2 ) “ ^ ("^l + TAtBvi) . 


vi IP + II V2 IP +rAt < Bvi,vi > +rAt < Bv2,V2 >) 


< ylui,U2 > — < ^U2,ui >= — (II vi Ip + II V2 iP +rAt < Bvi,vi > +rAt < Bv2,V2 >) 
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We now set for convenience N{vi,v2) =|| vi |p + || V 2 |P +rAt < Bvi,vi > +rAt < Bv 2 ,V 2 >■ 
We infer from the previous identities 


a = 1 — 


At 


N{vi,V2) 


{< Avi,vi > + < Av 2 ,V 2 >) and 6 = 


At 


N{vi,V2) 


< (^ - A^)vi,V2 > . 


The stability condition is 


■q = + h^ <1. 

After the usual simplifications, we find as a sufficient condition 


( 22 ) 


q < 1-2 


At^ 


At 

N{vi,V2) 


a (< Bvi,vi > + < Bv 2 , V 2 >) 


N{vi,V2 

At^ 


(< Bvi,Vi > + < Bv 2,V2 >y 


A-A' 


T ||2 


IP + II V2 IP)" 


' 4A(ui,U2p 

We now let Z =< Bvi,vi > + < Bv 2 ,V 2 > and Y =|| ui |p + || V 2 |p. The last inequality reads 


At {f3^ - 2ar)Z2 + 


A-A'^ IP ) < 2aYZ. 


1 


At this point, we set ^ ^ 5 I’ ^ ■ Note that we have \min{B) < ^ < 

^ II '^1 II + II ^2 II 

^max{B). After usual simplifications, the sufficient stability condition (22) writes now as 

At - 2aT)i + < 2a, 

where we have set 5 =|| A — A^ ||. We use the function <h(^) = (/3^ — 2ar)^ + ^<5^ ; is 

obviously regular on [Xmin{B),\max{B)], since Xmin{B) > 0, i? is assumed to be SPD. We 
deduce the following sufficient stability conditions 

• if <1?(^) < 0, G [Amm(-B), XmaxiB)], the scheme is unconditionnaly stable 

• if there exists ^ G [Xmin{B) , Xmax{B)\ such that <1>(^) > 0, then a sufficient stablility 
condition is 

-SfaP- WY 

\\min {B) i^max (B)] 


We conclude by studying the two cases. 

Unconditional stablity : To have <1>(^) < 0, G [XminiB), Xmax{B)], we must have ^{XminiB)) < 
0 and ^{Xmax{B)) < 0, that is 












We now remark that = (/3^-2Q;r)-^ < 0 r > ^ [^min{B), \max{B)\ 

. This is satisfied under the previous hypothesis, hence the first statement [z.] of the theorem. 


Conditional stability : The condition is Max(<l>(x)) > 0. We distinguish two cases 

,52 


<f> is monotone. 


- $'({) > 0 and > 0 ie r < 

o2 c: 

A sufficent stability condition is then r < ^ — -—r—^ 


and T < w - 




8aXmax{B) 


2 • 


0 < At < 


2a 8aXmax{Bf 
2a 


and 


^>(A max) 

- ^'{0 < 0 and HX^iniB)) > 0 that is ^ ^ < ll + SaxlniBf 


and 


0 < At < 


2a 


^Xmin{B))^ 
= 0 for ^ e]XminiB), Xm.ax{B)[. This means that 


^ 

2 q : 8aXmin{,B') 


/32 

<T < ^ - 


<52 


2cr 8aXraax(,B) 


2 • 


This implies that ^'{Xmin{B)) < 0 so Max(<l>(^)) is reached at Xmin{B) or at Xmax{B) and 
the stability condition is 


0 < At < 


2q: 


Max($(Amm(S)),^>(A max (B))) 

■ We point out that the symmetric part of matrix A is dominant for small values of 5 and 
that in this case the above stability result is to be compared with that of Proposition (2.1). 

r2 « 

Particularly inconditional stability condition is obtained for r > ^ + - —r-T—~ 4j, for 
inconditional preconditioners B as those presented above. 


4 Numerical Results 

4.1 The problem and the numerical schemes 

As stated in the introduction, we here aim at computing the steady state of the so-called 2D 
driven cavity problem in a rectangular domain XI. The steady state will be reached by a pseudo¬ 
temporal method, and for that purpose we consider the stream function-vorticity formulation 
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(w — V’) (see [ElEO] and the references therein): 


duj 

1 dcfidu dcpduj . 

(23) 


A'ljj = u, in n 

(24) 


0 

3 

II 

o' 

3 

(25) 


that we supplement with proper boundary conditions. We denote by Fj i = 1, ..,4 the sides of 
the unit square as follows: Fi is the lower horizontal side, F 3 is the upper horizontal side, F 2 
is the left vertical side, and F 4 is the right vertical side. 


U=g,V=0 



Figure 2: The lid driven cavity - Schematic localization of the mean vortex regions 

We distinguish two different driven flows, according to the choice of the boundary conditions 
on the velocity. More precisely we have 

• g{x) = 1 : Cavity A (lid driven cavity) 

• g{x) = (1 — (1 — 2x)^)^ : Cavity B (regularized lid driven cavity) 

We will also consider the steady linear part of this equation (the Stokes problem), whose the 
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solution will be chosen as the initial condition of (23) 


o' 

II 

3 

<1 

1 

1 

in U 

(26) 

w an= Alp ao 

on5U, 

(27) 

Alp = 00, 

in H 

(28) 

w 

II 

o 

on 

(29) 


The RSS scheme is based on two different finite differences discretization of differential operators 
at the same grid points: a second order finite difference scheme will be used for preconditiong 
while the fourth order compact scheme is implemented for the effective approximation to the 
solution. 

The boundary conditions on oj are derived by the discretization of on the boundaries. 
With the conditions on u and v we have 


uj{x, 0, t) 
uj(x, 1, t) 
u;(0,y,t) 
cv(l,y,t) 


Wi, 

dx^ 


(x,0,t) 

on Ti, 

(x,l,t) 

on Ts, 

(0,y,t) 

on r2, 

(i,y,i) 

on r4. 


So, since V’an = 0 and u = 


dip 

~dy' 



we obtain by using Taylor expansions 


<^i,0 




WOj 


WAT+lj 


_ 

2}? ’ 

^ - Qhgjih) 

21? ' 
_ - 8V^2,i , 

“ 2/i2 ’ 

_ —V^jv-1,7 + 

2/i2 


(30) 


Here g{x) denotes the boundary condition function for the horizontal velocity at the boundary 

Ta. 

The boundary conditions on ip are homogeneous Dirichlet BC. The operators are discretized by 
second order centered schemes on a uniform mesh composed by N points in each direction of 


the domain of step-size h = 


1 

N + 1 


The total number of unknowns is then 2N‘^. 


The boundary conditions on io are iteratively implemented according to the relations (30 30), 
making the finite differences scheme second order accurate. Using the following fourth order 
accurate extrapolation. 
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^1,0 

Wj,Ar+i 

^0,j 

WAT+lj 


1 / 8 1 \ 

- 3V'i,2 + gV’i.S - g^*,4 j , 

1 / 8 1 \ 25 

^ f 8V’i,Ar - 3V>i,Ar-l + gV’i,Ar-2 “ -'i/’i.Af-S j “ ^di^h), 

7:2 ( + gV'S.i - gV’ 4 J j , 


h? 

1 


1 


p- {^i>N,j - 3V’iV-l,j + g'0iV-2j - gV^Af-3,j ) , 


we complete the discretization. Now the implementation of the RSS scheme reads as 
• Convection-Diffusion problem: knowing compute solution of 


(^(fe+i) _ 1 (fc-i-i) 

Re ^ dy dx dx dy 


(31) 


At 


u 


(fe+i) = 


0 in D =]0, Ip 

on dQ 


(32) 


Poisson problem: knowing compute solution of 


j^(fc-i-i) _ ^^(fc+i) inD=]0,lp, 

^(fc-i-i) _ g onclD. 


(33) 


Algorithm 2 : RRS Navier-Stokes 


1 

2 

3 

4 

5 

6 


{io^,Tp^) given as solution of the Stokes problem (26) 
for A: = 0,1, • • • until convergence do 

Update the boundary terms using (31) 

Compute by solving. (32) with RSS ([^ 

Compute ^p'^~^^as solution of the Poisson equation (33) 
end for 


4.2 RSS Schemes for computing Steady States of the lid driven cavity 

We give now numerical results on the square cavity D =]0,1[^, we compare the numerical values 
of the steady state with those of the literature. Here g{x) = 1. 

The steady state is computed for ||^|| < e = 10“^. We report hereafter the vorticity and 
the stream function in figures Him 3jiici |6^ for Rg- — 100; R,g — 400; Rg — 1000; R,g — 3200 
respectively. They agree with those of the literature laiaiiaiis]; the localization of the extrema 
of ui and V' are reported on table 
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Omega 




Figure 3: Steady solution of NSE (23) - g' = 1 - r 


1 - iV = 127 - = 100 - At = 0.001 




Figure 4: Solution of NSE (23) - 5 = 1 - r 


1 - A = 127 - iie = 400 - At = 0.001 


Now we illustrate the influence of the stabilization parameter r on the convergence in time 
to the steady state. A large value of r allows to take a larra time step At but slows down 
the convergence in time. We have plotted in figures ([^ and (|8) the evolution in time of ||^||. 
We observe that, for r = 1 the RSS schemes (first and second order) behave similarly as the 
reference one (semi backward Euler’s); for r = 100, we see that RSS is slow downed while its 
extrapoled version has comparable dynamics to Euler’s. 

We now give numerical results for the rectangular cavity. They are presented in figures ([^ 
,([Tol), ( [IT| ) and (??) for Re = 100, Re = 400, Re = 1000, Re = 3200 respectiveley. They agree 
with those obtained by Goyon [13], see also the numerical values reported in Table 
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Omega 




Figure 6: Solution of NSE 


(23) 


g 


l - T = l - N = m - Re = 3200 - At = 0.0005 
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Re = 100 


Spatial accuracy 

Principal Vortex 

Ben Artzi et al |3] 

4th order compact scheme 

0. Goyon [T3] 
second order 

extrapoled RSS 
4th order 

grid 


97 X 97 

129 X 129 

127 X 127 

At 



0.005 

0.004 


intensity 


0.1033 

0.1026 


X 


0.6172 

0.6172 


y 


0.7343 

0.7422 


Re = 400 


Spatial accuracy 

Principal Vortex 

Ben Artzi et al j3] 

4th order compact scheme 

0. Goyon [T3] 
second order 

extrapoled RSS 
4th order 

grid 


97 X 97 

129 X 129 

127 X 127 

At 




0.017 


intensity 

0.1136 


0.1123 


X 

0.5521 


0.5625 


y 

0.6042 


0.6094 


Re = 1000 


Spatial accuracy 

Principal Vortex 

Ben Artzi et al |3] 

4th order compact scheme 

0. Goyon [13] 
second order 

extrapoled RSS 
4th order 

grid 


97 X 97 

129 X 129 

127 X 127 

At 



0.02 

0.01 


intensity 

0.1178 

0.1157 

0.1158 


X 

0.5312 

0.5312 

0.5391 


y 

0.5625 

0.5625 

0.5703 


Re = 3200 


Spatial accuracy 

Principal Vortex 

Ben Artzi et al |3| 

4th order compact scheme 

0. Goyon [13] 
second order 

extrapoled RSS 
4th order 

grid 


97 X 97 

129 X 129 

127 X 127 

At 



0.01 

0.006 


intensity 

0.1174 

0.1122 

0.1129 


X 

0.5208 

0.5234 

0.5156 


y 

0.5417 

0.5468 

0.5391 


Table 2: Extremal values of 'll: and their localization for NSE (23). Comparison with results of 
Croisille and Goyon, for different Reynolds numbers, g = 1 - t = 1 
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Figure 7: Convergence to NSE steady state (23)-r = l- A/' = 63--Re = 100 - At = 0.01 



Figure 8: Convergence to NSE steady state (23) 


r = 100 - iV = 63 - iie = 100 - At = 0.01 
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Figure 9: 


Figure 10: 
At = 0.001 



Solution of (23) in [0; 1] x [0; 2] - g 


1 - r = 1 - 255 X 511 - = 100 - At = 0.001 
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Figure 11: Solution of NSE ([2^ in [0; 1] x [0; 2] - 5 = 1 - r = 1 - 255 x 511 - Re = 1000 
At = 0.0005 
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Re = 100 




extrapoled RSS 

extrapoled RSS 

0. Goyon [13] 
second order 

C-H. Bruneau and 
C. Jouron |3] 

grid: 


63 X 127 

255 X 511 

65 X 129 

multi-grid 

At 


1 .10“^ 

1 .10“^ 

1 .10“^ 


e: 



RT^ 

Iir^ 


VS 


0.1040 

0.1034 

0.1035 

0.1033 


X 

0.6094 

0.6172 

0.6093 

0.6172 


y 

1.7344 

1.7344 

1.7343 

1.7344 

VI 

V' 

8 X 10"^ 

0.0003 

6.65 X 10-4 

0.783 X 10"^ 


X 

0.5469 

0.5820 

0.5468 

0.5391 


y 

0.5938 

0.5039 

0.5781 

0.5859 


Re = 400 




extrapoled RSS 

extrapoled RSS 

0. Goyon [T3| 
second order 

C-H. Bruneau and 
C. Jouron |3] 

grid: 


63 X 127 

255 X 511 

65 X 129 

multi-grid 

At 


1 .10“^ 

1 .10"^ 

1 .10"^ 


e: 



RT^ 

RT^ 


VS 

V' 

0.1131 

0.1120 

0.1097 

0.1124 


X 

0.5469 

0.5586 

0.5625 

0.5547 


y 

1.6094 

1.6094 

1.6094 

1.5938 

VI 

i’ 

0.009 

0.0059 

8.06 X 10“^ 

0.909 X 10“^ 


X 

0.4219 

0.5156 

0.4375 

0.4297 


y 

0.8438 

0.8750 

0.8593 

0.8125 


Re = 1000 




extrapoled RSS 

extrapoled RSS 

0. Goyon [T3| 
second order 

C-H. Bruneau and 
C. Jouron |3] 

grid: 


127 X 257 

255 X 511 

257 X 513 

multi-grid 

At 


5T0=^ 

5l0^^ 

5.10“^ 


e: 


RT^ 

RT^ 

RT^ 


VS 


0.1189 

0.1196 

0.1187 

0.1169 


X 

0.5312 

0.5312 

0.5313 

0.5273 


y 

1.5781 

1.5781 

1.5781 

1.5625 

VI 

V' 

0.0134 

0.0135 

1.32 X 10"^ 

0.0148 


X 

0.3438 

0.3398 

0.3359 

0.3516 


y 

0.8438 

0.8438 

0.8476 

0.7891 


Table 3: Some extremal values of uj and 
(upper and lower vortex) - g= l- T = l 


xjj for the steady state of NSE (23) on [0; 1] x [0; 2] 
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Figure 12; Solution of NSE (23) on [0; 1] x [0; 2] - 17 = 1 - r = 1 - 255 x 511 - Re = 3200 - 
At = 0.0005 


Here the vorticity uj and the streamfunction '0 at steady state for Re = 3200. The maximal 
values found for iIj are, for the maximum 0.0196 located in {x,y) = (0.4492,06914) and the 
minimum is —1.215, located in {x,y) = (0.5156,1.15664). 


4.3 On the role of r and of the extrapolation to the convergence in time to 
the steady state 

We now denote by Tc the time at which the convergence to the steady state is reached, say 

Tc = {k + 1 ) At with II -- II ^ We put the symbol * * * when the scheme is unstable 

and blows up numerically; ’NC’ when it is not convergent after A given large time (T=2000) 
and finally ’Inc. Stab’ when it is inconditionally stable. 

We first consider the RRS scheme with the second order laplacian matrix as a preconditioner. 
We show for low Reynolds numbers, that a large value of r allows to use large At and that the 
slower convergence to the steady state can be corrected by using the Richardson extrapolation. 


Of course, since the stabilization allows to take larger time steps, a important gain in CPU 
time can be obtained when computing a steady state. It can be estimated by considering the 
number of iteration in time NT: ^ for RSS and 3^ for RSS with extrapolation. For example, 
taking Re = 400 and n = 127, we find NT = 2501 for r = 1 and At = 0.014 and for r = 30 and 
At = 0.6, we have NT = 410 (RSS) and NT = 568 (RSS with extrapolation) ; hence a factor 
6.1 is reached for RSS and one of 4.4 for extrapolated RSS, see Table 
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T 

Extrapolation 

At 

no 

^^raax 

no 

Tc 

no 

At 

yes 

^tmax 

yes 

Tc 

yes 

r = 1 

0.01 

0.019 

15.62 

0.01 

0.019 

15.6 


0.019 


15.665 

0.019 


16.492 

r = 10 

0.02 

0.14 

16.86 

0.02 

0.11 

15.8 


0.06 


19.32 

0.06 


16.8 


0.07 


19.95 

0.07 


17.29 


0.1 


21.7 

0.1 


21.7 

r = 30 

0.3 

Inc. Stab. 

54.5 

0.3 

1.08 

32.1 


0.4 


79.6 

0.4 
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Table 4: RSS (left and RSS with Extrapolation (right) Re = 100, n = 63, e = 10 ^ 


T 

Extrapolation 

At 

no 

^^max 

no 

Tc 

no 

At 

yes 

^tmax 

yes 

Tc 

yes 

r = 1 

0.01 

0.014 

35.02 

0.01 

0.014 

35.03 


0.014 


35.014 

0.014 


35.042 

r = 20 

0.2 

0.28 

67.60 

0.21 

0.22 

57.6 

r = 30 

0.3 


105 

0.3 

0.35 

58.2 


0.35 


117.95 

0.35 


92.4 

r = 50 

0.3 

5.1 

139.2 

0.3 

1 

69 


0.4 


176.8 

0.4 


83.6 


0.5 


212.5 

0.5 


99 


0.6 


246.6 

0.6 


113.4 


Table 5: RSS (left) and RSS witht Extrapolation (right) Re = 400, n = 127, e = 10 ® 


We now consider larger Reynolds numbers and take into account the convective part of the 
equation in the construction of the sparse RSS preconditioner and apply nonlinear RSS scheme 
(|^ to the vorticity time marching, say 

A2 + diag{Dy^l)^^'^)Dx - diag{Dx'ip^^'^)Dii] 


ffc+l) _ Y 

At 


where A 2 is the second order laplacian matrix, diag{Dy'ijj^^^) (resp. diag{Dx'il>^^'^)) is the diagonal 

matrix with the discrete (second order accurate) approximation of (resp. ) at grid 

points as entries; Dx (resp. Dy) denote the (second order accurate) first derivative matrix in 
X (resp. in y) on the cartesian grid. Finally is the high order compact scheme 


discretisation of — + 


dcj) doj 


_ _ d(f) du} 

Re ~dy ~5x J)x IhJ' 

As Shown on Table and Table NLRSS outperform RSS for the computation of Steady 
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Method 

RSS 

RSS 

RSS 

RSS 

RSS 

RSS 

NLRSS 

NLRSS 

NLRSS 

T 

At 


Tc 

At 

^tmax 

Tc 

At 

^tmax 

Tc 

Extrapolation 

no 

no 

no 

yes 

yes 

yes 

yes 

yes 

yes 

r = 1 

0.005 

0.005 

56.21 

0.005 

0.01 

56.81 





0.01 



0.01 

56.79 


0.01 

0.02 

56.86 


0.02 


*** 

0.02 

*** 


0.02 


56.96 

T = 30 

0.05 

0.04 

NC 

0.05 

0.08 

47.95 

0.05 

0.7 

65.05 


0.1 


*** 

0.1 


*** 

0.1 


62.5 


0.7 


*** 

0.7 


*** 

0.7 


321.3 


Table 6: RSS (left) RSS with Extrapolation (center) and extrapolated NLRSS (right) Re = 1000, 
n = 127, e = lO'^ 


r 

Extrapolation 

At 

no 

^^max 

no 

Tc 

no 

At 

no 

^^max 

no 

Tc 

no 

r = 10 

0.1 

0.6 

223.9 

0.1 

0.006 



Table 7; NLRSS with (left) and RSS (right) Re = 3200, n = 127, e = 10 ® 


States for higher Reynolds numbers, for Re = 1000 and a fortiori for Re = 3200. It allows to 
use large times steps while RSS is unstable for such At. 

5 Concluding remarks 

We have studied RRS-like scheme (and their implementations) and pointed out their advantages 
for the numerical solution of parabolic problems when using high order compact schemes in 
finite differences for the space disctretization. In particular, the possibility of using fast solvers 
attached to a standard second order discretization, speeds up the resolution while bringing an 
enhanced stability. We also pointed out the role of the approximation of tB to A in the dynamics 
of the convergence to a steady state: a too strong stablization slows down the convergence in time 
while enhacing the stability of the scheme, the application of Richardson extrapolation allows 
to recover a dynamics close to the one of the classical schemes. The robustness of the schemes is 
illustrated with the solution of 2D NSE equations. The RSS approach is very versatile and allows 
adaptations of a large number of techniques of numerical analysis of ODEs. Many developments 
remain to consider, such as applying factorization updatings on the preconditioners or deriving 
and applying multilevel general (or Block) RSS schemes for the solution of other large scale 
parabolic problems. 
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